clear all

sysdir set PLUS ado_rdrobust760/plus
which rdrobust
which rdbwselect2

version 14.2
set type double
set matsize 8000
global current_path = c(pwd)
global path0 = "../firmleveltest"

cap mkdir "$path0"
cap mkdir "$path0/testresult_020119"
cap mkdir "$path0/graph_sp1x2"
cap mkdir "$path0/graph_combined"
cap mkdir "$path0/graph_combined_pdf"
cap mkdir "$path0/graph_combined_gph"

****************************************************************************************
use "../data/biddata020119.dta", clear
merge 1:1 rid using "../data/rid_address_license_firmID_020119.dta"
order firmID
drop _merge

/*** ad_firm no_license はoriginalのもの、ad_firm_newは付け替えをやったもの*/
replace ad_firm = ad_firm_new if ad_firm=="" & ad_firm_new ~= ""
replace no_license = no_license_new if no_license == . & no_license_new < .
drop *_new

set seed 1
set type double

keep if FLAG == 0
/* test for each firm and for 3 fiscal years*/
*drop samples including outlines*
bys id: egen min_tech = min(sc_tech)
drop if min_tech < 100
bys id: egen max_tech = max(sc_tech)
drop if max_tech > 250
/* drop auctions in which rank1_p bidder's q is 100*/
bys id: egen ddd = max (( rank1_p == 1 ) * (sc_tech == 100))
*drop if ddd == 1
drop ddd
/* drop bidders whose sc_tech == 100 in the loop */
*drop if sc_tech == 100

	
*** temporally drop unnecessary samples (before2012) ****
keep if fy_con > 2014 & fy_con < 2018
tab fy_con
********************************************************


local d = string( date("$S_DATE", "DMY"),"%td_CY-N-D")
di "`d'"

forvalues z = 1(1)2{
	cap drop *dep`z'*
	cap drop rank_i`z'*
	/* rank1_x > `z'を作る*/
	gen dep`z's = ldelta`z'p_s
	gen dep`z'p = ldelta`z's_p
	gen indep`z's = ldelta`z's_s 
	gen indep`z'p = ldelta`z'p_p
	replace dep`z's = -ldelta`z's_s if rank1_s == `z'
	replace dep`z'p = -ldelta`z'p_p if rank1_p == `z'
}


keep firmID nm_firm region *dep* no_license ad_firm rid id rank* fy_con *_date 
compress 

********* Tie auctions (score and price) excluded from RDD ******************
cap drop tie*
foreach k in "s" "p"{
	forvalues z = 1(1)2{
		gen aaa = 1 if indep`z'`k' == 0
		bys id: egen tie`z'`k' = max(aaa)
		drop aaa
	}
	replace tie2`k' = 1 if tie1`k' == 1
}

order nm_firm ad_firm no_license fy_con
format %50s nm_* ad_*
sort nm_firm ad_firm no_license

/******* to be deleted *****************/
cap drop groupregionfirm
egen groupregionfirm = group(firmID)
replace groupregionfirm = groupregionfirm * 10 + region
***************************************************** 

compress

graph drop _all
cap drop regionfirmID
gen regionfirmID = firmID * 10 + region
egen regionfirmTG = tag(regionfirmID)

********* For deleting obs which are too close to zero ****************
forvalues x = 1(1)2{
	cap drop rank_i`x'*
	bys rank1_s: egen rank_i`x's = rank(indep`x's) if rank1_s == `x' & indep`x's ~= 0
	bys rank1_p: egen rank_i`x'p = rank(-indep`x'p) if rank1_p == `x' & indep`x'p ~= 0
}	

cap drop ftt_*
forvalues t = 1(1)2{
	gen ftt_flip`t' = .
	gen ftt_flip`t'1 = .
	gen ftt_N`t' = .
	gen ftt_N`t'1 = .
	foreach k in "s" "p"{
		gen ftt_tau`t'`k' = .
		gen ftt_se`t'`k' = .
		gen ftt_ts`t'`k' = .
		gen ftt_rank1obsdrop`t'`k' = .
		foreach d in "l" "r"{
			gen ftt_bw`d'`t'`k' = .
			gen ftt_N`d'`t'`k' = .
			gen ftt_h`d'`t'`k' = .
		}
	}
}

compress
save $path0/sample020119.dta, replace

use $path0/sample020119.dta, clear
*pause on
*log using loglog2.log, replace
sum regionfirmID
local counter = 1
forvalues y = 1(1)`r(max)'{
	local no_partic = 0
	local quad_count = 0 /* counts # successful RD for each firm*/
	qui count if regionfirmID == `y'
	local no_partic = r(N)
	if `no_partic' > 4{
		di "[`counter']" _continue
		local counter = `counter' + 1
		cap drop theauction
		qui bys id: egen theauction = max( regionfirmID == `y' )
		cap replace thebidder = .
		if _rc ~= 0{
			qui gen thebidder = 1 if regionfirmID == `y' & regionfirmTG == 1
		}
		else{
			qui replace thebidder = 1 if regionfirmID == `y' & regionfirmTG == 1
		}
		preserve
			qui keep if theauction == 1|rank1_s <= 2 | rank1_p <= 2		
			sort thebidder
			local nmf = nm_firm[1]
			local fid = regionfirmID[1]
			local fad = ad_firm[1]
			local license = no_license[1]
			local region = region[1]
			
			/*if ustrregexm("`nmf'", "(?:三共建設（株）|岩崎電気（株）|（株）安井建設|北紘建設（株）|三井住友建設（株）|（株）安藤・間|砺波工業（株）|開成工業（株）|小林建設（株）|高砂建設（株）)") == 1 ///
			 & ustrregexm("`fad'",  "(?:北海道士別市東２北４－１５|東京都中央区日本橋馬喰町１－４－１６|長野県佐久市取出町４５０－５|北海道伊達市元町７１－２１|東京都中央区佃２－１－６|東京都港区赤坂６－１－２０|富山県砺波市中央町３－２１|香川県綾歌郡宇多津町浜一番丁１|広島県庄原市東本町４‐１‐１|三重県四日市市高砂町８－２９)") == 1 ///
			 {
			 	di "***********************"
			 }
			 else{
			 	restore
			 	continue
			 }
			*/
			************************************************************
			/* need to retry 
			三共建設（株）	北海道士別市東２北４－１５
			岩崎電気（株）	東京都中央区日本橋馬喰町１－４－１６
			（株）安井建設	長野県佐久市取出町４５０－５
			北紘建設（株）	北海道伊達市元町７１－２１
			三井住友建設（株）	東京都中央区佃２－１－６
			（株）安藤・間	東京都港区赤坂６－１－２０
			砺波工業（株）	富山県砺波市中央町３－２１
			開成工業（株）	香川県綾歌郡宇多津町浜一番丁１
			小林建設（株）	広島県庄原市東本町４‐１‐１
			高砂建設（株）	三重県四日市市高砂町８－２９
			*/
			************************************************************
			
				

			forvalues z = 1(1)2{
				if `z' == 2{
					qui keep if theauction == 1 | rank1_s == 2 | rank1_p == 2
				}
				qui count if dep`z's < 0 & indep`z's < 0 & theauction == 1
				local flip`z' = `r(N)'
				qui count if dep`z's < 0 & indep`z's < 0 & rank1_s == `z' + 1 ///
					& theauction == 1
				local flip`z'1 = `r(N)'
				qui count if indep`z's < 0 & theauction == 1
				local all`z' = `r(N)'
				qui count if indep`z's < 0 & rank1_s == `z' + 1 & theauction == 1
				local all`z'1 = `r(N)'
				foreach k in "s" "p" {
					local est`z'`k' = .
					local se`z'`k' = .
					local ts`z'`k' = .
					local rank1obsdrop`z'`k' = .  
					foreach d in "l" "r"{
						local bw`d'`z'`k' = .
						local h`d'`z'`k' = .
						local N`d'`z'`k' = .
					}
					
				
				************** whether x variable is score or price ****************************
					if "`k'" == "s"{
						local leftright = "l"
						/* rank1_s == 1のobsのうち、最も0に近いもので落とす数*/
						*local list_dropobs_trial = "1800 1200 800 200 1 0"
						local list_dropobs_trial = "3000 2900 2800 2700 2600 2500 2400 2300 2200 2100 2000 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100 90 80 70 60 50 40 30 20 10"
					}
					else if "`k'" == "p"{
						local leftright = "r"
						/* rank1_p == 1のobsのうち、最も0に近いもので落とす数*/
						*local list_dropobs_trial = "2000 1400 1000 888 500 2 1"
						local list_dropobs_trial = "3000 2900 2800 2700 2600 2500 2400 2300 2200 2100 2000 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100 90 80 70 60 50 40 30 20 10"

					}
				********************************************************************************
					*pause	
						/* Do not include auctions in which (tie`z'`k' == 1) highest scores/prices are tie*/
							/* scoreが同率の1位と2位の者の価格bidは通常違うので、一位をどちらの業者にするかで2位以下の者すべての価格乖離率dep`z'sが変わり、RDの結果が変わるのを防ぐため
							なお、set seed 1では解決にならないのは、サンプルが変わって、行の番号が変わると同率の1位と2位の者のどちらをrank1_s = 1にするかは変化するから*/
					local goloop = 1	
					foreach dr of numlist `list_dropobs_trial' {
						if `goloop' == 1{
							*di "z:`z';k:`k'; dr: `dr' @@@@@@@@@ "
							cap /*noisily*/ rdbwselect2 dep`z'`k' indep`z'`k' `leftright' if tie`z'`k' ~= 1 ///
								 & (theauction == 1 | (rank1_`k' == `z' & rank_i`z'`k' > `dr') ) /// /*note: rank_i1s <= `dr' excludes `dr' obs closest to zero for rd computation*/
									, kernel(tri) vce(cluster id) bwselect(msetwo)
							cap /*noisily*/ rdrobust dep`z'`k' indep`z'`k' if tie`z'`k' ~= 1 &  ///
								 (theauction == 1 | (rank1_`k' == `z' & rank_i`z'`k' > `dr') )  ///
									, kernel(tri) vce(cluster id) ///
								h(`e(h_msetwo_l)' `e(h_msetwo_r)') b(`e(b_msetwo_l)' `e(b_msetwo_r)') 
							if _rc == 0{
								local quad_count = `quad_count' + 1
								local goloop = 0
								local rank1obsdrop`z'`k' = `dr' /* this local records how many obs of rank1 closest to zero dropped*/

								if "`k'" == "s" {
									local est`z'`k' = string(`e(tau_bc_l)', "%6.3f" )
									mat B = e(V_rb_l)
									local se`z'`k' = string(sqrt(B[1,1]), "%6.3f" )
									local ts`z'`k' = string(`e(tau_bc_l)' / sqrt(B[1,1]) , "%6.3f" )
								}
								else if "`k'" == "p"{
									local est`z'`k' = string(`e(tau_bc_r)', "%6.3f" )
									mat B = e(V_rb_r)
									local se`z'`k' = string(sqrt(B[1,1]), "%6.3f" )
									local ts`z'`k' = string(`e(tau_bc_r)' / sqrt(B[1,1]) , "%6.3f" )
								}							
								foreach d in "l" "r" {
									local bw`d'`z'`k' = e(b_`d')
									local h`d'`z'`k' = e(h_`d')
									local N`d'`z'`k' = e(N_`d')
									*di "##### bandwidth2: `bw`d'`z'`k'' ######"
								}
							}
							else{
								*di as error "indep`z'`k' in (`fid') (obs.`no_partic') has an invertibility problem!"
							}
						}
					}
				}
				if `z' == 1 {
					/* ldelta1s_s*/
					di as result %-10s "(`fid')" %49s "`nmf';" /*%6.0f `fid' "-" `region'*/ ///
						"; Ns`z'=" %6.0g `Nl`z's' /// "; Fl=" %6.3f `flip`z'' / `all`z'' /// 
						"; est:" %6.3f `est1s' "(" %5.3f `se1s' "); t=" %6.3f `ts1s' "; Drobs:`rank1obsdrop`z's'"
					/* ldelta1p_p*/
						di as text " LID:" %12.0g `license' %42s "`fad'" ///
						"; Np`z'=" %6.0g `Nr`z'p' /// "; Fl=" %6.3f `flip`z'' / `all`z'' /// 
						"; est:" %6.3f `est1p' "(" %5.3f `se1p' "); t=" %6.3f `ts1p' "; Drobs:`rank1obsdrop`z'p'"
					*}
				}
				if `z' == 2 {
					/* ldelta1s_s*/
					di " # participations:" %5.0g `no_partic' ";" _continue
					di as result %36s " " %6.0f `xxx' /*"-" `region'*/ ///
						" Ns`z'=" %6.0g `Nl`z's' /// "; Fl=" %6.3f `flip`z'' / `all`z'' /// 
						"; est:" %6.3f `est2s' "(" %5.3f `se2s' "); t=" %6.3f `ts2s' "; Drobs:`rank1obsdrop`z's'"
					/* ldelta1p_p*/
						di as text %58s " " %10.0g " " %6.0f `xxx' /*"-" `region'*/ /// ///
						"; Np`z'=" %6.0g `Nr`z'p' /// "; Fl=" %6.3f `flip`z'' / `all`z'' /// 
						"; est:" %6.3f `est2p' "(" %5.3f `se2p' "); t=" %6.3f `ts2p' "; Drobs:`rank1obsdrop`z'p'"
					*}
				}
			}
		
			qui keep if thebidder == 1
			qui keep nm_firm ad_firm no_license firmID region* ftt_*

			forvalues z = 1(1)2{
				foreach k in "s" "p" {
					qui replace ftt_tau`z'`k' = `est`z'`k'' /* if regionfirmID == `y'*/
					qui replace ftt_se`z'`k' = `se`z'`k'' /* if regionfirmID == `y'*/
					qui replace ftt_ts`z'`k' = `ts`z'`k'' /* if regionfirmID == `y'*/
					qui replace ftt_rank1obsdrop`z'`k' = `rank1obsdrop`z'`k'' /* if regionfirmID == `y'*/
					foreach d in "l" "r" {
						qui replace ftt_bw`d'`z'`k' = `bw`d'`z'`k'' /* if regionfirmID == `y'*/
						qui replace ftt_N`d'`z'`k' = `N`d'`z'`k'' /* if regionfirmID == `y'*/
						qui replace ftt_h`d'`z'`k' = `h`d'`z'`k'' /* if regionfirmID == `y'*/
					}
				}
				qui replace ftt_N`z' = `all`z'' /* if regionfirmID == `y'*/
				qui replace ftt_N`z'1 = `all`z'1' /* if regionfirmID == `y'*/
				qui replace ftt_flip`z' = `flip`z'' /* if regionfirmID == `y'*/
				qui replace ftt_flip`z'1 = `flip`z'1' /* if regionfirmID == `y'*/
			}
			qui gen partic = `no_partic'
			cap /*noisily*/ save ///
			"$path0/testresult_020119/`y'_`nmf'_XN`no_partic'_Obs1s(`Nl1s'_`Nr1s')_Obs1p(`Nl1p'_`Nr1p')_pll.dta" ///
			, replace
		restore
	}
	
}


/* combining result files*/
clear all

cap erase $path0/testresults_020119.dta
save $path0/testresults_020119.dta, empty

local list : dir "$path0/testresult_020119" files "*.dta"
foreach x in `list' {
	*sleep 50
	di as text "Appending `x' file"
	append using "$path0/testresult_020119/`x'"
}
save $path0/testresults_020119.dta, replace


/* Draw graphs for each firm listed on resultsyyyy-mm-dd.dta*/
cd $current_path

use $path0/testresults_020119.dta, clear
*************************************************
drop regionfirmTG
merge 1:m regionfirmID using $path0/sample020119.dta
format %30s nm* ad*
**************
sort ftt_flip1
gen ftt_fr1 = ftt_flip1 / ftt_N1
cap drop tsrank
cap drop mints1
cap drop neg_ftt_
gen neg_ftt_ts1s = -ftt_ts1s
egen mints1 = rowmin(neg_ftt_ts1s ftt_ts1p)
egen aaa = rank(mints1) if regionfirmTG == 1, t
egen bbb = rank(regionfirmID)
gen abab = aaa + 1/bbb
egen bcbc = rank(abab) if regionfirmTG == 1
bys regionfirmID: egen tsrank = max(bcbc)
drop aaa bbb abab bcbc 
compress

cap drop _merge
save "$path0/test_sample_with_results.dta", replace


/* it takes a couple of hours to draw all graphs*/
use "$path0/test_sample_with_results.dta", clear
/* creating firm list with excel format*/
order nm_firm ad_firm no_license firmID region mints1 tsrank ftt_tau1s ftt_ts1s ftt_tau1p ftt_ts1p ftt_Nl1s ftt_Nr1p 
sort mints1

merge 1:1 rid using "../data/biddata020119.dta"
drop _merge
merge 1:1 rid using "../data/rid_address_license_firmID_020119.dta"
drop _merge
keep if FLAG == 0
bys id: egen min_tech = min(sc_tech)
bys id: egen max_tech = max(sc_tech)
drop if min_tech < 100 & max_tech > 100 
drop if max_tech > 250					/*drop price-only auctions*/

*************************************************************************
keep mints1 ftt* nm_firm ad_firm *firmID rank1_p rank1_s *drop* ///
	ldelta* fy_con no_license region mints1 tsrank id tie* regionfirmTG
	order *_s *_p nm_firm *drop* 
count if (( ftt_ts1s > 1.64 & ftt_ts1s < .)) & regionfirmTG == 1
count if ((ftt_ts1p < -1.64 & ftt_ts1p < . ))  & regionfirmTG == 1
gen firm_positive = ((ftt_ts1p < -1.64 & ftt_ts1p < . )|( ftt_ts1s > 1.64 & ftt_ts1s < .)) 
count if firm_positive == 1& regionfirmTG == 1
local num_positive_firms = r(N)


local _firm_tsranking = 1
local graph_counter = 1
local number_of_panels = 12
local ranking_init_graph = `_firm_tsranking'

local _graphs = ""
cap graph drop _all
di "num_positive_firms: `num_positive_firms'"
qui do "modules/draw_graph_firmleveltest1x2.do"
local continue = 1
while `continue' == 1 {
	*draw_graph_2x2 0 `_firm_tsranking' "" "$path0/graph/" "3" "`graph_counter'"
	*cap noisily graph drop firm_`graph_counter'_graph

	draw_graph_1x2 0 `_firm_tsranking' "" "$path0/Graph_sp1x2/" "3" "`graph_counter'"
	qui graph dir
	local _graphs_dir = r(list)
	if ustrregexm("`_graphs_dir'", "firm_`graph_counter'_graph"){
		if `_firm_tsranking' > `num_positive_firms'{
			cap graph drop firm_`graph_counter'_graph
			twoway (function y=0, lwidth(none) ytitle("") xtitle("") xscale(lc(white)) ///
			yscale(lc(white)) ylabel(,labcolor(white) nogrid noticks) ///
			xlabel(,labcolor(white) nogrid noticks)), name(firm_`graph_counter'_graph) xsize(10) ysize(5)
		}
		local _graphs = "`_graphs' firm_`graph_counter'_graph"
		local graph_counter = `graph_counter' + 1
	}
	if `graph_counter' == `number_of_panels' + 1{
		local ranking_last_graph = min(`_firm_tsranking',`num_positive_firms')
		local combined_graph_title = "firms_`ranking_init_graph'-`ranking_last_graph'"
		graph combine `_graphs', cols(2) rows(6) xsize(7) ysize(10) scale(.8) imargin(vsmall) ///
		title(`combined_graph_title',size(small)) graphregion(color(white)) holes(`_holes')
		*qui graph export "$path0/graph_combined/`combined_graph_title'.png", width(2750) replace
		qui graph save "$path0/graph_combined_gph/`combined_graph_title'.gph", replace
		*qui graph export "$path0/graph_combined_pdf/`combined_graph_title'.pdf", as(pdf) replace
		local graph_counter = 1
		graph drop _all
		local _graphs = ""
		local ranking_init_graph = `_firm_tsranking' + 1
		if `_firm_tsranking' >= `num_positive_firms'{
			local continue = 0
		}
	}
	local _firm_tsranking = `_firm_tsranking' + 1
}	

global path_gph ../firmleveltest/graph_combined_gph
global path_png ../firmleveltest/graph_combined_pdf

forvalues x = 1(12)1141{
	local x12 = `x' + 11
	if `x' == 1141{
		local x12 = `x' + 2
	}
	local filename = "firms_`x'-`x12'"
	if fileexists("$path_png/`filename'.pdf"){
		di "`filename'.pdf exists"
	}
	else{
		qui graph use "$path_gph/`filename'"
		cap noisily graph export "$path_png/`filename'.pdf", as(pdf)
		sleep 200
		graph drop _all
	}
}
